We’ll carry on looking at spatial data in this lab. You should have the following files from the previous lab, but it’s worth checking to make sure you have everything:
Base R has no structure for spatial data, so you will need to install the following packages (you should have some of these from previous modules):
library(ggplot2)
library(terra)
library(RColorBrewer)
library(sf)
library(viridis)
Previously, we were working with vector spatial data. These are geometries composed of points defined by their coordinates. An alternative form of spatial data is known as a raster. This is gridded data. It takes the form of a rectangle composed of squares of equal size, which are sometimes called ‘cells’ or ‘pixels’. Each cell stores some kind of value.
This simplifies the geometry, which can be specified by two pieces of
information: the spatial extent of the raster and the resolution of the
cells. Here we create a blank raster with 10 rows and 10 columns, with a
resolution of 10x10 using the rast() function. We then
assign random values to each cell:
r <- rast(nrow = 10, ncol = 10,
xmin = 0, xmax = 100,
ymin = 0, ymax = 100)
r[] <- runif(n = 100)
r
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source(s) : memory
## name : lyr.1
## min value : 0.004815464
## max value : 0.972897647
rast objects can be plotted using the base
plot() command:
plot(r)
The terra package offers a wide array of functions
for dealing with gridded data, including the ability to read from many
widely used file formats, like remote sensing images (e.g. GeoTiffs),
NetCDF, and HDF formats. We will use it here to work with a Landsat 8
scene collected on June 14, 2017. The subset covers the area between
Concord and Stockton, in California, USA. The files are contained in the
zip file rs.zip. Download this (if you haven’t already) and
unzip it in your data folder.
We will also need the shapefile of California places (ca_places.zip), so download and unzip this as well. We’ll read this in before starting:
ca_places <- st_read("./data/ca_places/ca_places.shp")
## Reading layer `ca_places' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/ca_places/ca_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1618 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.2695 ymin: 32.53432 xmax: -114.229 ymax: 41.99317
## Geodetic CRS: WGS 84
To read in gridded data, use the rast() function. This
will read in the Near Infrared (NIR) channel (B5)
b5 <- rast("./data/rs/LC08_044034_20170614_B5.tif")
b5
## class : SpatRaster
## dimensions : 1245, 1497, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source : LC08_044034_20170614_B5.tif
## name : LC08_044034_20170614_B5
## min value : 0.0008457669
## max value : 1.0124315023
When we print the rast object created, the second line
of the output lists the dimensions of the data. Note that here, this has
1245 rows, 1497 columns and 1 layer. This also shows the resolution
(30x30 m), the extent, the CRS and the a brief summary of the data.
We can write rast objects back to file using
writeRaster() (I’ll bet you never thought it would be
called that). You can write out to any format supported by GDAL. Here we
write out to a TIFF format. You can see the full list of available
formats for reading and writing by running the
writeFormats() function in the console. We’ll use this
again after having worked with the data.
writeRaster(b5,
filename = "./b5.tif",
overwrite = TRUE)
As with the sf objects, we can check the coordinate
reference system of the file we just read. This does not print very
well, but if you look at the end you’ll see a reference to the EPSG code
(32610).
crs(b5)
## [1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 10N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-123,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n BBOX[0,-126,84,-120]],\n ID[\"EPSG\",32610]]"
If the CRS is not set, you can set it using the crs()
function and an EPSG code. For example, the following code would set the
CRS to WGS84 (don’t run this as the CRS is already
defined)
crs(b5) <- "EPSG:4326"
Again, this should not be used to change the CRS, only set
it. Note that crs is for rast objects,
st_crs for vectors.
You can transform the CRS for a raster layer using
project(). This can again use a EPSG type code.
b5_wgs84 <- project(b5, "EPSG:4326")
crs(b5_wgs84)
## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
We’ll keep the Landsat data in its UTM projection (10N), and reproject the CA place data to match:
ca_places <- st_transform(ca_places, 32610)
You can make a simple plot using the plot()
function:
plot(b5, main = "Landsat 8 (B2)")
plot(st_geometry(ca_places), add = TRUE)
The function cellStats() can be used to calculate most
summary statistics for a raster layer. So to get the mean global
temperature (and standard deviation):
global(b5, mean)
## mean
## LC08_044034_20170614_B5 0.2693644
global(b5, sd)
## sd
## LC08_044034_20170614_B5 0.111683
If we want to use only a subset of the original raster layer, the
function crop() will extract only the cells in a given
region. This can be defined using another raster object or Spatial*
object, or by defining an extent object:
# extent method
my_ext <- ext(c(xmin = 612500,
xmax = 617500,
ymin = 4196000,
ymax = 4201000))
b5_sub <- crop(b5, my_ext)
plot(b5_sub)
We can also use an sf object to crop the data. Here,
we’ll extract the polygon for Bethel Island from the
ca_places object, and use this to crop the raster:
bethel <- ca_places %>%
dplyr::filter(NAME == "Bethel Island")
b5_sub <- crop(b5, bethel)
plot(b5_sub)
plot(st_geometry(bethel), add = TRUE)
Note that crop subsets the original raster to the extent
of Canada’s borders, rather than to the borders themselves. This is
because rasters are always rectangular. You can ‘hide’ the
values of raster cells outside of a polygon by using the
mask function. The raster has to be rectangular, so this
does not remove the cells outside the polygon. Rather, it sets their
value to NA.
b5_sub <- mask(b5_sub, mask = bethel)
plot(b5_sub)
plot(st_geometry(bethel), add = TRUE)
Values can be extracted from individual locations (or sets of
locations) using extract(). This can take a set of
coordinates in matrix form, or use a Spatial* object. To get the
reflectance value at 615000 E, 4199000 N:
extract(b5, cbind(615000,4199000))
## LC08_044034_20170614_B5
## 1 0.2121801
You can also extract for multiple locations. Let’s generate a set of random points in Bethel Island, and then sample the reflectance value for each of these.
random_pnts <- st_sample(bethel, size = 20)
extract(b5, st_coordinates(random_pnts))
## LC08_044034_20170614_B5
## 1 0.10249011
## 2 0.27481058
## 3 0.35472512
## 4 0.20094655
## 5 0.06612195
## 6 0.09268785
## 7 0.21324277
## 8 0.32954717
## 9 0.36669603
## 10 0.46638858
## 11 0.41353872
## 12 0.31954971
## 13 0.39768595
## 14 0.34171325
## 15 0.46441513
## 16 0.48172089
## 17 0.36231536
## 18 0.39916062
## 19 0.34785050
## 20 0.28463453
You can also extract values within a polygon, by replacing the point
coordinates with a sf polygon:
b5_bethel <- extract(b5, bethel)
head(b5_bethel)
## ID LC08_044034_20170614_B5
## 1 1 0.10251180
## 2 1 0.08084705
## 3 1 0.05189565
## 4 1 0.17219034
## 5 1 0.29502234
## 6 1 0.31579795
By default, this returns the value of all pixels within the polygon.
By adding the fun= argument, you can easily calculate
summary statistics:
extract(b5, bethel, fun = 'median')
## ID LC08_044034_20170614_B5
## 1 1 0.3779947
Note that if the sf object has multiple polygons, it
will return the summary statistic for each one. Let’s now extract the
values for a different place (Oakley), and then we can compare them. We
do this in a couple of steps. First get the polygon for Oakley, then
extract using the two polygons combined.
oakley <- ca_places %>%
dplyr::filter(NAME == "Oakley")
b5_bethel_oakley <- extract(b5, rbind(bethel, oakley))
names(b5_bethel_oakley) <- c("ID", "B5")
We’ll visualize the difference with ggplot2, which shows higher NIR reflectance values for Bethel, likely indicating higher vegetation cover.
ggplot(b5_bethel_oakley, aes(x = B5, fill = as.factor(ID))) +
geom_histogram(alpha = 0.7, position = 'identity')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The b2 raster represents a single band from the Landsat
8 scene. More usefully, we can load several bands and combine them into
a single object. Here, we’ll load the blue (B2), green (B3), red (B4),
and infrared (B5) bands :
# Blue
b2 <- rast('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- rast('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- rast('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- rast('data/rs/LC08_044034_20170614_B5.tif')
Now we’ll create a raster stack with all 4 of these:
s <- c(b5, b4, b3, b2)
s
## class : SpatRaster
## dimensions : 1245, 1497, 4 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## sources : LC08_044034_20170614_B5.tif
## LC08_044034_20170614_B4.tif
## LC08_044034_20170614_B3.tif
## LC08_044034_20170614_B2.tif
## names : LC08_04~0614_B5, LC08_04~0614_B4, LC08_04~0614_B3, LC08_04~0614_B2
## min values : 0.0008457669, 0.02084067, 0.04259216, 0.0748399
## max values : 1.0124315023, 0.78617686, 0.69246972, 0.7177562
The metadata now shows that this contains 4 layers. You can also make the stack directly by passing a list of file names:
filenames <- paste0('./data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- rast(filenames)
landsat
## class : SpatRaster
## dimensions : 1245, 1497, 11 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## sources : LC08_044034_20170614_B1.tif
## LC08_044034_20170614_B2.tif
## LC08_044034_20170614_B3.tif
## ... and 8 more source(s)
## names : LC08_~14_B1, LC08_~14_B2, LC08_~14_B3, LC08_~14_B4, LC08_~14_B5, LC08_~14_B6, ...
## min values : 0.09641791, 0.0748399, 0.04259216, 0.02084067, 0.0008457669, -0.007872183, ...
## max values : 0.73462820, 0.7177562, 0.69246972, 0.78617686, 1.0124315023, 1.043204546, ...
This now contains all bands representing reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.
If we now use the extract() function with this object,
it will return values for all bands in the stack:
extract(landsat, cbind(615000,4199000))
## LC08_044034_20170614_B1 LC08_044034_20170614_B2 LC08_044034_20170614_B3
## 1 0.1691325 0.1547111 0.1395956
## LC08_044034_20170614_B4 LC08_044034_20170614_B5 LC08_044034_20170614_B6
## 1 0.1370149 0.2121801 0.1594604
## LC08_044034_20170614_B7 LC08_044034_20170614_B8 LC08_044034_20170614_B9
## 1 0.1526508 0.09177701 0.001236123
## LC08_044034_20170614_B10 LC08_044034_20170614_B11
## 1 311.6736 308.671
par(mfrow = c(2,2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))
plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))
The values in each layer range from 0 to 1, and the same scale has been used for each band, showing clearly the difference in reflectance for this set of wavelengths. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter. In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark
The bands can be combined to form composite images. Here, we’ll use
the red, green and blue bands make a true color image. This uses the
concatenation function (c()), and the order of the bands is
important (R, then G, then B):
landsatRGB <- c(b4, b3, b2)
plotRGB(landsatRGB, stretch = "lin")
We can also make a false color composite with the NIR, red and green bands, where bright reds indicate vegetation cover:
Another popular image visualization method in remote sensing is known “false color” image in which NIR, red, and green bands are combined. This representation is popular as it makes it easy to see the vegetation (in red).
landsatFCC <- c(b5, b4, b3)
plotRGB(landsatFCC, stretch="lin")
# Raster algebra
terra makes it easy to carry out simple raster algebraic operations. Eahc band or layer is treated a 2D array which makes it possible to add, subtract, multiply, divide, etc. As a simple example here, we can calculate the NDVI for the Landsat scene as
\[ NDVI = (NIR - R) / (NIR + R) \]
NIR is band 5 and red is band 4:
ndvi <- (b5 - b4) / (b5 + b4)
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")
We can also make a quick histogram of values to look for any outliers:
hist(ndvi, main = "NDVI")
## Warning: [hist] a sample of54% of the cells was used
We can then use the results to carry out some simple classification.
clamp function masks
all value otuside of a range (here below 0.4)veg = clamp(ndvi, lower=0.4, values=FALSE)
plot(veg)
crops = ndvi > 0.25 & ndvi < 0.3
plot(crops)
And water bodies (NDVI < 0):
water = ndvi < 0
plot(water)
In the previous lab, we looked at making maps in base R (with the
plot() function and an sf object) as well as
using ggplot2. More recently, there has been something
of an explosion in mapping packages, including interactive maps. We’ll
look here at two of these: tmap and
leaflet. Other useful ones include
mapsf highcharter,
mapview and mapdeck.
library(sf)
## Census tracts for Salt Lake County with population density
tracts <- st_read("./data/slc_tract/slc_tract_2015.shp")
## Reading layer `slc_tract_2015' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/slc_tract/slc_tract_2015.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 212 features and 17 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -112.2602 ymin: 40.41417 xmax: -111.5532 ymax: 40.92185
## Geodetic CRS: NAD83
## Salt Lake light rail tracks
lightrail <- st_read("./data/LightRail_UTA/LightRail_UTA.shp")
## Reading layer `LightRail_UTA' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRail_UTA/LightRail_UTA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 5 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 413241.3 ymin: 4486347 xmax: 429455.1 ymax: 4515254
## Projected CRS: NAD83 / UTM zone 12N
## Salt Lake light rail stations
stations <- st_read("./data/LightRailStations_UTA/LightRailStations_UTA.shp")
## Reading layer `LightRailStations_UTA' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRailStations_UTA/LightRailStations_UTA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 57 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 413245.8 ymin: 4486448 xmax: 429389.5 ymax: 4515197
## Projected CRS: NAD83 / UTM zone 12N
## Landsat images for California
# Blue
b2 <- rast('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- rast('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- rast('data/rs/LC08_044034_20170614_B4.tif')
## Places for California
ca_places <- st_read("./data/ca_places/ca_places.shp")
## Reading layer `ca_places' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/ca_places/ca_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1618 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.2695 ymin: 32.53432 xmax: -114.229 ymax: 41.99317
## Geodetic CRS: WGS 84
## W North America Climate
wna_climate <- read.csv("./data/WNAclimate.csv")
## Natural Earth country shapefiles
countries <- st_read("./data/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")
## Reading layer `ne_50m_admin_0_countries' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 241 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS: WGS 84
## Natural Earth country rivers
tmap makes thematic maps. It works in a very similar
way to ggplot2, by building a series of layers and map
geometries and elements. In general, we start by using
tm_shape() to identify the spatial object to be used, and
add geometries are added, including filled polygons, borders and
symbols. Finally, we can add legends, scale bars, etc. Note that you
will need to install this package if you haven’t done so already
(install.packages("tmap")).
We’ll first make a map with a mixture of vector layers, including census tracts, light rail tracks and stations. Start by loading these
library(sf)
## Census tracts for Salt Lake County with population density
tracts <- st_read("./data/slc_tract/slc_tract_2015.shp")
## Reading layer `slc_tract_2015' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/slc_tract/slc_tract_2015.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 212 features and 17 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -112.2602 ymin: 40.41417 xmax: -111.5532 ymax: 40.92185
## Geodetic CRS: NAD83
## Salt Lake light rail tracks
lightrail <- st_read("./data/LightRail_UTA/LightRail_UTA.shp")
## Reading layer `LightRail_UTA' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRail_UTA/LightRail_UTA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 5 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 413241.3 ymin: 4486347 xmax: 429455.1 ymax: 4515254
## Projected CRS: NAD83 / UTM zone 12N
## Salt Lake light rail stations
stations <- st_read("./data/LightRailStations_UTA/LightRailStations_UTA.shp")
## Reading layer `LightRailStations_UTA' from data source
## `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRailStations_UTA/LightRailStations_UTA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 57 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 413245.8 ymin: 4486448 xmax: 429389.5 ymax: 4515197
## Projected CRS: NAD83 / UTM zone 12N
If we check the CRS, you’ll see that tracts has a
different one, so will need to be reprojected:
st_crs(tracts)$epsg
## [1] 4269
st_crs(lightrail)$epsg
## [1] 26912
tracts <- st_transform(tracts, st_crs(lightrail))
Now let’s start mapping things out with tmap. First,
let’s make a simple map showing the polygon outlines using
tm_borders(). This is the standard tmap
method: first define the object that you want to plot with
tm_shape, then define what you want to plot.
library(tmap)
tm_shape(tracts) +
tm_borders()
The function tm_fill() will then fill these using one of
the variables in the tract data set. We’ll use the population density
(density). Note that this automatically adds a legend
within the frame of the figure:
tm_shape(tracts) +
tm_fill("density") +
tm_borders()
The color scale can be changed by setting the palette
argument in tm_fill(). This includes the ColorBrewer scales
described in the previous module, and the algorithm used to define the
intervals. For example, to use the ‘Greens’ palette with percentile
breaks. We also change the color of the tract boundaries to
lightgray:
tm_shape(tracts) +
tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray")
You can reverse the color scale by prepending a -
(-Greens). If you want to see the full set of palettes that
you can use, install the tmaptools package and run the
following code (you may also need to install
shinyjs):
tmaptools::palette_explorer()
Let’s now add another layer. We’ll overlay the light rail track on
this map. As this is a different object in R (lightrail),
we first need to use tm_shape to indicate that we are using
it, then add tm_lines() to show the routes. We’ll make the
lines a little thicker with lwd and change the type to
dashed.
tm_shape(tracts) +
tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2, lty = 'dashed', col = "darkorange")
Note that if you use a varible name for the color, it will thematically map the lines (I’ve removed the polygon fill to make this clear):
tm_shape(tracts) +
#tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 4, col = "ROUTE")
We can now add the stations (using
tm_shape() to select the
correct object to plot):
tm_shape(tracts) +
tm_fill("density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23)
Now let’s add some details to the map. We’ll add a graticule before
the polygons are plotted, a compass and a scalebar. We’ll also use the
tm_layout function to add a title and move the legend to a
new position. (This function has a lot of options for changing your map
- it’s well worth checking the help page.)
tm_shape(tracts) +
tm_graticules(col = "lightgray") +
tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23) +
tm_compass(position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "top")) +
tm_layout(main.title = "Salt Lake County Light Rail",
legend.outside = TRUE,
legend.outside.position = c("left"))
The map can be clipped to smaller regions by using the
bbox argument in the first tm_shape()
function. In this code, we’ll extract a set of census tracts that
correspond to downtown Salt Lake. We’ll then use the the bounding box of
these to crop the map (I’ve removed the fill). We’ll also add the
station name using tm_text():
tracts_sub = tracts %>%
dplyr::filter(TRACTCE %in% c(102500, 102600, 114000))
tm_shape(tracts, bbox = st_bbox(tracts_sub)) +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23) +
tm_text("STATIONNAM", ymod = -1, bg.color = "white", size = 0.8)
Finally, tmap allows you to make interactive maps
fairly simply with the tmap_mode() function. Setting this
to view will make an interactive map, to plot
will make a static map (like the ones we have made so far). Here we’ll
remake the full map as an interactive map:
## Set interactive
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(tracts) +
tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile", id = "TRACTCE") +
tm_borders("lightgray") +
tm_shape(lightrail) +
tm_lines(lwd = 2) +
tm_shape(stations) +
tm_dots(size = 0.25, shape = 23)
## Symbol shapes other than circles or icons are not supported in view mode.
We’ll reset back to static maps for the rest of this exercise:
tmap_mode("plot")
## tmap mode set to plotting
In this example, we’ll map out the climate data from western North
America. We’ll go a bit faster in this example and not explain every
step. First read in the data and create an sf object:
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
wna_climate <- read.csv("./data/WNAclimate.csv")
wna_climate <- st_as_sf(wna_climate,
coords = c("LONDD", "LATDD"),
crs = 4326)
Make a quick thematic map of average july temperature:
tm_shape(wna_climate) +
tm_symbols(col = "Jul_Tmp")
Next, we’ll download country outlines and river centerlines from Natural Earth using the
rnaturalearth package (you will need to install this).
This will download layers from the Natural Earth website - there are a
variety of these available, including country and state boundaries,
physical features (rivers, lakes, etc) and cultural features (place
names, roads, etc). There are three levels of resolution: 1:10m, 1:50m
and 1:110m. To download a layer, you need to specify the resolution
(scale) and the name of the layer. The full set can be
found here
library(rnaturalearth)
countries50 <- ne_download(scale = 50, type = "admin_0_countries")
## Reading layer `ne_50m_admin_0_countries' from data source
## `/private/var/folders/ql/nw995vq50pq3dlrxhk7mm4_40000gq/T/RtmpiMUnCe/ne_50m_admin_0_countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 242 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS: WGS 84
rivers50 <- ne_download(scale = 50, type = "rivers_lake_centerlines", category = "physical")
## Reading layer `ne_50m_rivers_lake_centerlines' from data source
## `/private/var/folders/ql/nw995vq50pq3dlrxhk7mm4_40000gq/T/RtmpiMUnCe/ne_50m_rivers_lake_centerlines.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 478 features and 36 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -165.2439 ymin: -50.24014 xmax: 176.3258 ymax: 73.3349
## Geodetic CRS: WGS 84
Now we can put this all together as a series of layers. We’ll make the color palette continuous and use a red-blue palette.
tm_shape(countries50, bbox = st_bbox(wna_climate)) +
tm_borders() +
tm_shape(rivers50) +
tm_lines("lightblue", lwd = 2) +
tm_shape(wna_climate) +
tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
style = "cont", title.col = "degC")
## Warning: The shape rivers50 contains empty units.
Finally, we’ll make two maps (july temperature and annual
precipitation), and use tmap_arrange() to make a two panel
plot:
m1 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
tm_borders() +
tm_shape(rivers50) +
tm_lines("lightblue", lwd = 2) +
tm_shape(wna_climate) +
tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
style = "cont", title.col = "degC")
m2 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
tm_borders() +
tm_shape(rivers50) +
tm_lines("lightblue", lwd = 2) +
tm_shape(wna_climate) +
tm_symbols(col = "annp", palette = "BuPu", alpha = 0.75,
style = "cont", title.col = "mm/yr")
tmap_arrange(m1, m2)
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
## Warning: The shape rivers50 contains empty units.
In the final example, we’ll use the Landsat data with
tmap. The main function to map raster data is
tm_raster(). We’ll use it here to show the ndvi layer made
earlier.
tm_shape(ndvi) +
tm_raster()
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
We’ll update this a little by clamping the NDVI values at a lower limit of 0, changing the color palette and making it continuous. We’ll also add the California places added as an overlay.
tm_shape(clamp(ndvi, lower = 0)) +
tm_raster(palette = "Greens", style = "cont", title = "NDVI") +
tm_shape(ca_places) +
tm_borders(lwd = 2) +
tm_layout(main.title = "Landsat 8 (2017/06/14)",
legend.position = c("left", "top"),
legend.bg.color = "white", legend.bg.alpha = 0.7)
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)
A second function (tm_rgb()) allows the creation of
color composites. We’ll use the red, green and blue bands to make the
true color composite. The maximum reflectance value in these files is
1.0, so we need to set this for scaling.
tm_shape(landsat) +
tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)
## Warning: Raster values found that are outside the range [0, 1]
When we made this figure earlier, it was a lot brighter. This is
because the previous function ‘stretched’ the range of values for each
layer. We can mimic that here by first using the stretch()
function from terra, then repeating the plot:
landsat_stretch = stretch(landsat, minv = 0, maxv = 1,
minq = 0.01, maxq = 0.99)
tm_shape(landsat_stretch) +
tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)
The tmap mode function allows you to quickly make an
interacitve map. The resulting map is based on the well-known Leaflet
javascript library, which offers a much greater range of flexibility in
making these maps. Fortunately, there is an R library
(leaflet) that allows access to the Leaflet API
allowing you to build quite complex maps from R. We’ll remake the map of
Salt Lake County and the light rail routes using this now. First install
the library (install.packages("leaflet")), and load it (and
a helper library for working with html formats):
library(leaflet)
library(htmltools)
leaflet operates in a similar way to
ggplot2 and tmap in building a map by
layers. Layers are added using the pipe operator that we encountered in
an earlier lab (%>%). There are three basic ingredients
that are required:
leaflet() functionHere’s a really simple example, where we use the default tiles from OpenStreetMap and drop a marker on the University of Utah:
m <- leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m # Print the map
You can also try changing the Provider - the source of the background map. For example, this will use ESRI map (the map is not shown here but will appear in RStudio when you run this):
m <- leaflet() %>%
addProviderTiles("Esri.WorldStreetMap") %>%
addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m
ESRI imagery:
m <- leaflet() %>%
addProviderTiles("Esri.WorldImagery") %>% # Add default OpenStreetMap map tiles
addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m
You can find a full set of these providers here: https://leaflet-extras.github.io/leaflet-providers/preview/
We’ll now add the census tracts and light rail sf
objects. First, these need to be reprojected to WGS84 long-lat (EPSG
code 4326) to work with Leaflet.
tracts_ll <- st_transform(tracts, crs = 4326)
stations_ll <- st_transform(stations, crs = 4326)
lightrail_ll <- st_transform(lightrail, crs = 4326)
Now we can add these to a dark basemap. First the tracts with
addPolygons, then the routes (addPolylines)
and stations (addMarkers). Note that the stations include a
label taken from the station name column in the
sf object. This wil appear when you hover your cursor over
a marker.
leaflet() %>%
# add a dark basemap
addProviderTiles("CartoDB.DarkMatter") %>%
# add the polygons of the clusters
addPolygons(
data = tracts_ll,
color = "#E2E2E2",
opacity = 1, # set the opacity of the outline
weight = 1, # set the stroke width in pixels
fillOpacity = 0.2 # set the fill opacity
) %>%
addPolylines(
data = lightrail_ll
) %>%
addMarkers(
data = stations_ll,
label = ~htmlEscape(STATIONNAM)
)
Finally, we’ll add some extra details. First, we’ll set colors for each station by the light rail line it lies on. This is a little complicated - below we:
sf objectunique(stations_ll$LINENAME)
## [1] "University" "N/S TRAX" "Med Center"
## [4] "Intermodal Hub" "Airport" "Draper"
## [7] "Mid-Jordan" "West Valley" "Sugar House Streetcar"
line_pal <- RColorBrewer::brewer.pal(9, "Set1")
line_no <- as.numeric(as.factor(stations_ll$LINENAME))
stations_ll$LINECOL <- line_pal[line_no]
Next we’ll make a popup window for each station. This will appear
when the marker is clicked (as opposed to the label which appears when a
cursor hovers over the marker). This has to be formatted using html.
Here, we use R’s paste() function to stick together html
tags with the station name, address, etc. The <b>
tags make the labels bold, and the <br> tags add a
line break. You can do a lot more here, for example, including images or
URLs in each popup.
station_popup = paste0(
"<b>Station: </b>",
stations_ll$STATIONNAM,
"<br>",
"<b>Line Name: </b>",
stations_ll$LINENAME,
"<br>",
"<b>Park n Ride: </b>",
stations_ll$PARKNRIDE,
"<br>",
"<b>Address: </b>",
stations_ll$ADDRESS
)
With all this in place, we can rebuild our map. The main changes here are that we use a circle marker as this has an option to change color, and include the popup information in the marker.
leaflet() %>%
# add a dark basemap
addProviderTiles("Esri.WorldStreetMap") %>%
# add the polygons of the clusters
addPolygons(
data = tracts_ll,
color = "#E2E2E2",
# set the opacity of the outline
opacity = 1,
# set the stroke width in pixels
weight = 1,
# set the fill opacity
fillOpacity = 0.2
) %>%
addPolylines(
data = lightrail_ll
) %>%
addCircleMarkers(
data = stations_ll,
color = ~LINECOL,
label = ~htmlEscape(STATIONNAM),
popup = station_popup
)